rm(list = ls())

## Load libraries

library(ggplot2, warn.conflicts = FALSE)
library(data.table, warn.conflicts = FALSE)
library(plm, warn.conflicts = FALSE)
suppressPackageStartupMessages(library(tidyverse, warn.conflicts = FALSE))

## Load data

load("../generated_data/bes_w1_w21_pacer_w1_w2_merged.Rdata")

setkey(bes_long, id, wave)

########################
## Reshape data for analysis

## Outcome variables

outcome_vars <- c("redistSelf", "ptvCon", "ptvLab", "ptvLD", "taxSpendSelf" ,"trustWestminster" ,"deficitReduce", "riskTaking", "approveUKGovt", "govtHandleCostLive", "govtHandleEcon", "govtHandleImmig", "govtHandleNHS", "govtHandleEduc", "govtHandleLevelCrime", "govtCorpSupport", "ideoBHPSPprivateEnterprise", "ideoBHPSSstateOwnership", "ideoBHPSJjobForAll", "ideoNoNeedUnions", "ideoFairShareWealth", "ideoOneLawRich", "reasonForUnemployment", "immigrantsWelfareState", "govtHandouts", "benefitsNotDeserved", "riskPoverty", "riskUnemployment", "lockdownApproval")

bes_long <- tibble(bes_long)

out <- bes_long %>%
  filter((id_in_any_pacer == TRUE|id_in_bes_w20 == TRUE) &
           id_in_wave == TRUE) %>%
  filter(wave < 21) %>%
  mutate(income_guarantee_treatment = ifelse(is.na(income_guarantee_binary),FALSE, income_guarantee_binary),
         income_guarantee_covid_treatment = ifelse(is.na(income_guarantee_covid_binary),FALSE, income_guarantee_covid_binary),
         income_guarantee_covid_job_treatment = ifelse(is.na(income_guarantee_covid_job_binary),FALSE, income_guarantee_covid_job_binary),
         income_guarantee_covid_self_treatment = ifelse(is.na(income_guarantee_covid_self_binary),FALSE, income_guarantee_covid_self_binary),
         income_guarantee_other_treatment = ifelse(is.na(income_guarantee_other_binary),FALSE, income_guarantee_other_binary),
         workingStatus_imputed_ordinal = factor(workingStatus,
                                                  levels = rev(c("Working full time (30 or more hours per week)", 
                                                                 "Working part time (8-29 hours a week)", 
                                                                 "Working part time (less than 8 hours a week)", 
                                                                 "Unemployed and looking for work"))),
         workingStatus_imputed_ordinal_2 = factor(recode(workingStatus,
                                                         "Not in paid work for any other reason" = "Unemployed and looking for work"),
                                                  levels = rev(c("Working full time (30 or more hours per week)", 
                                                                 "Working part time (8-29 hours a week)", 
                                                                 "Working part time (less than 8 hours a week)", 
                                                                 "Unemployed and looking for work"))),
         vote19 = recode_factor(p_past_vote_2019,
                                Conservative = "Conservative",
                                Labour = "Labour",
                                .default = "Other"),
         income_guarantee_furlugh_ws = case_when(is.na(workingStatus) ~ NA,
                                                 workingStatus == "Furloughed" ~ TRUE,
                                                 workingStatus != "Furloughed" ~ FALSE),
         polAttention_numeric = as.numeric(factor(polAttention, levels = c("Pay no attention",1:9,"Pay a great deal of attention")))) %>%
  rename(cross_sectional_weight = cross_sectional_weight_) %>%
  group_by(id) %>%
  mutate(lost_employment = ifelse(wave < 19.5 | length(wave[wave<19.5]) == 0, FALSE, 
                                  as.numeric(workingStatus_imputed_ordinal)[wave == max(wave[wave < 19.5])] > 
                                    as.numeric(workingStatus_imputed_ordinal)),
         lost_employment = ifelse(is.na(lost_employment), FALSE, lost_employment),
         lost_employment_2 = ifelse(wave < 19.5 | length(wave[wave<19.5]) == 0, FALSE, 
                                    as.numeric(workingStatus_imputed_ordinal_2)[wave == max(wave[wave < 19.5])] > 
                                      as.numeric(workingStatus_imputed_ordinal_2)),
         lost_employment_2 = ifelse(is.na(lost_employment_2), FALSE, lost_employment_2),
         
         # Pre-crisis attitudes
         ## Mean prior attitudes
         redistSelf_prior = mean(as.numeric(redistSelf[wave <= 19]), na.rm = TRUE),
         taxSpendSelf_prior = mean(as.numeric(taxSpendSelf[wave <= 19]), na.rm = TRUE),
         
         ## Mean prior attitudes - BES14 onwards
         redistSelf_prior14 = mean(as.numeric(redistSelf[wave <= 19 & wave >= 14]), na.rm = TRUE),
         taxSpendSelf_prior14 = mean(as.numeric(taxSpendSelf[wave <= 19 & wave >= 14]), na.rm = TRUE),
         
         ## Recent prior attitudes
         redistSelf_prior_recent = ifelse(sum(wave<=19 & !is.na(redistSelf)) != 0,
                                          as.numeric(redistSelf[wave == max(wave[wave<=19 & !is.na(redistSelf)])]),
                                          NA),
         
         taxSpendSelf_prior_recent = ifelse(sum(wave<=19 & !is.na(taxSpendSelf)) != 0,
                                          as.numeric(taxSpendSelf[wave == max(wave[wave<=19 & !is.na(taxSpendSelf)])]),
                                          NA),
         
         ## Don't know attitudes
         redistSelf_dk_prior = any(redistSelf_dk[wave <= 19], na.rm = TRUE),
         taxSpendSelf_dk_prior = any(taxSpendSelf_dk[wave <= 19], na.rm = TRUE),
         
         ### Remove any observations for which we don't have any prior attitudes
         none_redistSelf_prior = all(is.na(redistSelf_dk[wave <= 19])),
         none_taxSpendSelf_prior = all(is.na(taxSpendSelf_dk[wave <= 19])),
         redistSelf_dk_prior = ifelse(none_redistSelf_prior, NA, redistSelf_dk_prior),
         taxSpendSelf_dk_prior = ifelse(none_taxSpendSelf_prior, NA, taxSpendSelf_dk_prior),
         
         ## Political attention
         polAttention_prior = mean(as.numeric(polAttention_numeric[wave <= 19]), na.rm = TRUE)
         
         ) %>%
  ungroup() 

# Specify prior attitude high and low groups

## Averaging across all pre-crisis waves

median_redistSelf_prior <- out %>% group_by(id) %>% summarise(x = unique(redistSelf_prior)) %>% pull(var = x) %>% median(na.rm = TRUE)
median_taxSpendSelf_prior <- out %>% group_by(id) %>% summarise(x = unique(taxSpendSelf_prior)) %>% pull(var = x) %>% median(na.rm = TRUE)
median_polAttention_prior <- out %>% group_by(id) %>% summarise(x = unique(polAttention_prior)) %>% pull(var = x) %>% median(na.rm = TRUE)

## Averaging across all pre-crisis waves from BES14 onwards
median_redistSelf_prior14 <- out %>% group_by(id) %>% summarise(x = unique(redistSelf_prior14)) %>% pull(var = x) %>% median(na.rm = TRUE)
median_taxSpendSelf_prior14 <- out %>% group_by(id) %>% summarise(x = unique(taxSpendSelf_prior14)) %>% pull(var = x) %>% median(na.rm = TRUE)

## Most recently expressed prior attitude
median_redistSelf_prior_recent <- out %>% group_by(id) %>% summarise(x = unique(redistSelf_prior_recent)) %>% pull(var = x) %>% median(na.rm = TRUE)
median_taxSpendSelf_prior_recent <- out %>% group_by(id) %>% summarise(x = unique(taxSpendSelf_prior_recent)) %>% pull(var = x) %>% median(na.rm = TRUE)


out <- out %>%
  mutate(polAttention_prior_group = factor(ifelse(polAttention_prior >= median_polAttention_prior, "High", "Low"), levels = c("Low", "High")),
         redistSelf_prior_group = case_when(redistSelf_prior <= 5 ~ "Right",
                                            redistSelf_prior >=5 & redistSelf_prior <=7  ~ "Mid",
                                            redistSelf_prior >= 7 ~ "Left"),
         redistSelf_prior_group14 = case_when(redistSelf_prior14 <= 5 ~ "Right",
                                              redistSelf_prior14 >=5 & redistSelf_prior14 <=7  ~ "Mid",
                                              redistSelf_prior14 >= 7 ~ "Left"),
         redistSelf_prior_group_recent = case_when(redistSelf_prior_recent <= 5 ~ "Right",
                                                   redistSelf_prior_recent >=5 & redistSelf_prior_recent <=7  ~ "Mid",
                                                   redistSelf_prior_recent >= 7 ~ "Left"),
         
         taxSpendSelf_prior_group = case_when(taxSpendSelf_prior <= 5 ~ "Right",
                                              taxSpendSelf_prior >=5 & taxSpendSelf_prior <=7  ~ "Mid",
                                              taxSpendSelf_prior >= 7 ~ "Left"),
         
         taxSpendSelf_prior_group14 = case_when(taxSpendSelf_prior14 <= 5 ~ "Right",
                                                taxSpendSelf_prior14 >=5 & taxSpendSelf_prior14 <=7  ~ "Mid",
                                                taxSpendSelf_prior14 >= 7 ~ "Left"),
         
         taxSpendSelf_prior_group_recent = case_when(taxSpendSelf_prior_recent <= 5 ~ "Right",
                                                     taxSpendSelf_prior_recent >=5 & taxSpendSelf_prior_recent <=7  ~ "Mid",
                                                     taxSpendSelf_prior_recent >= 7 ~ "Left")
  ) 


############################
## Outputs for draft

## Function for implementing plm on arbitrary outcomes

plm_1 <- function(y, my_formula = NULL){
  #cat(paste0(y," \n"))
  
  out$Y <- scale(as.numeric(out[[y]]))
  
  plm(as.formula(my_formula),
      effect = "twoway",
      model = "within",
      index = c("id","wave"),
      data = out,
      weights = out$cross_sectional_weight) 
}


## Define models

## - Baseline spec = collapse UC and furlough

outcome_vars <- c("redistSelf", "taxSpendSelf", "ideoOneLawRich", "ideoBHPSJjobForAll", "reasonForUnemployment", "riskPoverty", "riskUnemployment", "govtHandouts")

baseline_formula <- "Y ~ lost_employment + income_guarantee_treatment"

baseline_out_income_guarantee <- lapply(outcome_vars, function(x) plm_1(y = x, my_formula = baseline_formula)) 

names(baseline_out_income_guarantee) <- outcome_vars

### - Baseline spec = separating UC, furlough, self-employment

outcome_vars <- c("redistSelf", "taxSpendSelf", "ideoOneLawRich", "ideoBHPSJjobForAll", "reasonForUnemployment", "govtHandouts")

baseline_formula <- "Y ~ lost_employment + income_guarantee_covid_job_treatment + income_guarantee_covid_self_treatment + income_guarantee_other_treatment"

baseline_out_income_guarantee_separate <- lapply(outcome_vars, function(x) plm_1(y = x, my_formula = baseline_formula)) 

names(baseline_out_income_guarantee_separate) <- outcome_vars


## - party interaction spec = separating UC, furlough, self-employment

outcome_vars <- c("redistSelf", "taxSpendSelf", "ideoOneLawRich", "ideoBHPSJjobForAll", "reasonForUnemployment", "govtHandouts")

baseline_formula <- "Y ~ lost_employment + income_guarantee_covid_job_treatment*vote19 + income_guarantee_covid_self_treatment*vote19 + income_guarantee_other_treatment*vote19 + vote19*as.factor(wave)"

party_out_income_guarantee_separate <- lapply(outcome_vars, function(x) plm_1(y = x, my_formula = baseline_formula)) 

names(party_out_income_guarantee_separate) <- outcome_vars


## - prior position interaction spec = separating UC, furlough, self-employment

baseline_formula <- "Y ~ lost_employment*redistSelf_prior_group + income_guarantee_covid_job_treatment*redistSelf_prior_group + income_guarantee_covid_self_treatment*redistSelf_prior_group + income_guarantee_other_treatment*redistSelf_prior_group + redistSelf_prior_group * as.factor(wave)"

prior_out_income_guarantee_separate_redist <- plm_1(y = "redistSelf", my_formula = baseline_formula)

baseline_formula <- "Y ~ lost_employment*taxSpendSelf_prior_group + income_guarantee_covid_job_treatment*taxSpendSelf_prior_group + income_guarantee_covid_self_treatment*taxSpendSelf_prior_group + income_guarantee_other_treatment* taxSpendSelf_prior_group+ taxSpendSelf_prior_group * as.factor(wave)"

prior_out_income_guarantee_separate_taxspend <- plm_1(y = "taxSpendSelf", my_formula = baseline_formula)

prior_out_income_guarantee_separate <- list(redistSelf = prior_out_income_guarantee_separate_redist,
                                            taxSpendSelf = prior_out_income_guarantee_separate_taxspend)

## - prior position interaction spec, w14 onwards = separating UC, furlough, self-employment

baseline_formula <- "Y ~ lost_employment*redistSelf_prior_group14 + income_guarantee_covid_job_treatment*redistSelf_prior_group14 + income_guarantee_covid_self_treatment*redistSelf_prior_group14 + income_guarantee_other_treatment*redistSelf_prior_group14 + redistSelf_prior_group14* as.factor(wave)"

prior_out_income_guarantee_separate_redist14 <- plm_1(y = "redistSelf", my_formula = baseline_formula)

baseline_formula <- "Y ~ lost_employment*taxSpendSelf_prior_group14 + income_guarantee_covid_job_treatment*taxSpendSelf_prior_group14 + income_guarantee_covid_self_treatment*taxSpendSelf_prior_group14 + income_guarantee_other_treatment*taxSpendSelf_prior_group14 + taxSpendSelf_prior_group14 * as.factor(wave)"

prior_out_income_guarantee_separate_taxspend14 <- plm_1(y = "taxSpendSelf", my_formula = baseline_formula)

prior_out_income_guarantee_separate14 <- list(redistSelf = prior_out_income_guarantee_separate_redist14,
                                            taxSpendSelf = prior_out_income_guarantee_separate_taxspend14)

## - prior position interaction spec, most recent prior position = separating UC, furlough, self-employment

baseline_formula <- "Y ~ lost_employment*redistSelf_prior_group_recent + income_guarantee_covid_job_treatment*redistSelf_prior_group_recent + income_guarantee_covid_self_treatment*redistSelf_prior_group_recent + income_guarantee_other_treatment*redistSelf_prior_group_recent + redistSelf_prior_group_recent * as.factor(wave)"

prior_out_income_guarantee_separate_redist_recent <- plm_1(y = "redistSelf", my_formula = baseline_formula)

baseline_formula <- "Y ~ lost_employment*taxSpendSelf_prior_group_recent + income_guarantee_covid_job_treatment*taxSpendSelf_prior_group_recent + income_guarantee_covid_self_treatment*taxSpendSelf_prior_group_recent + income_guarantee_other_treatment*taxSpendSelf_prior_group_recent + taxSpendSelf_prior_group_recent * as.factor(wave)"

prior_out_income_guarantee_separate_taxspend_recent <- plm_1(y = "taxSpendSelf", my_formula = baseline_formula)

prior_out_income_guarantee_separate_recent <- list(redistSelf = prior_out_income_guarantee_separate_redist_recent,
                                              taxSpendSelf = prior_out_income_guarantee_separate_taxspend_recent)


## - don't know interaction spec = separating UC, furlough, self-employment

baseline_formula <- "Y ~ lost_employment*redistSelf_dk_prior + income_guarantee_covid_job_treatment*redistSelf_dk_prior + income_guarantee_covid_self_treatment*redistSelf_dk_prior + income_guarantee_other_treatment*redistSelf_dk_prior + redistSelf_dk_prior * as.factor(wave)"

dk_out_income_guarantee_separate_redist <- plm_1(y = "redistSelf", my_formula = baseline_formula)

baseline_formula <- "Y ~ lost_employment*taxSpendSelf_dk_prior + income_guarantee_covid_job_treatment*taxSpendSelf_dk_prior + income_guarantee_covid_self_treatment*taxSpendSelf_dk_prior + income_guarantee_other_treatment*taxSpendSelf_dk_prior + taxSpendSelf_dk_prior * as.factor(wave)"

dk_out_income_guarantee_separate_taxspend <- plm_1(y = "taxSpendSelf", my_formula = baseline_formula)

dk_out_income_guarantee_separate <- list(redistSelf = dk_out_income_guarantee_separate_redist,
                                        taxSpendSelf = dk_out_income_guarantee_separate_taxspend)


## - political attention interaction spec = separating UC, furlough, self-employment

baseline_formula <- "Y ~ lost_employment*polAttention_prior_group + income_guarantee_covid_job_treatment*polAttention_prior_group + income_guarantee_covid_self_treatment*polAttention_prior_group + income_guarantee_other_treatment*polAttention_prior_group + polAttention_prior_group * as.factor(wave)"

attention_out_income_guarantee_separate_redist <- plm_1(y = "redistSelf", my_formula = baseline_formula)

attention_out_income_guarantee_separate_taxspend <- plm_1(y = "taxSpendSelf", my_formula = baseline_formula)

attention_out_income_guarantee_separate <- list(redistSelf = attention_out_income_guarantee_separate_redist,
                                                taxSpendSelf = attention_out_income_guarantee_separate_taxspend)


## - prior government support interaction spec = separating UC, furlough, self-employment

baseline_formula <- "Y ~ lost_employment*gov_support_past + income_guarantee_covid_job_treatment*gov_support_past + income_guarantee_covid_self_treatment*gov_support_past + income_guarantee_other_treatment*gov_support_past + gov_support_past * as.factor(wave)"

prior_gov_support_out_income_guarantee_separate_redist <- plm_1(y = "redistSelf", my_formula = baseline_formula)

prior_gov_support_out_income_guarantee_separate_taxspend <- plm_1(y = "taxSpendSelf", my_formula = baseline_formula)

prior_gov_support_out_income_guarantee_separate <- list(redistSelf = prior_gov_support_out_income_guarantee_separate_redist,
                                                taxSpendSelf = prior_gov_support_out_income_guarantee_separate_taxspend)

## - interacting lost_employment with benefit receipt

outcome_vars <- c("redistSelf", "taxSpendSelf", "ideoBHPSJjobForAll", "riskUnemployment", "riskPoverty")

baseline_formula <- "Y ~ lost_employment * income_guarantee_other_treatment + income_guarantee_covid_job_treatment + income_guarantee_covid_self_treatment "

baseline_out_income_guarantee_separate_lost_employment_interaction <- lapply(outcome_vars, function(x) plm_1(y = x, my_formula = baseline_formula)) 

names(baseline_out_income_guarantee_separate_lost_employment_interaction) <- outcome_vars


## Save all model outputs

save(out, 
     baseline_out_income_guarantee, 
     baseline_out_income_guarantee_separate, 
     party_out_income_guarantee_separate, 
     prior_out_income_guarantee_separate,
     prior_out_income_guarantee_separate14,
     prior_out_income_guarantee_separate_recent,
     dk_out_income_guarantee_separate,
     attention_out_income_guarantee_separate,
     prior_gov_support_out_income_guarantee_separate,
     baseline_out_income_guarantee_separate_lost_employment_interaction,
     file = "../generated_data/fe_models_new.Rdata")
